// P1Diffractor.cpp: implementation of the P1Diffractor class.
//
//////////////////////////////////////////////////////////////////////

#include "Pbar1Diffractor.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

Pbar1Diffractor::Pbar1Diffractor(char *DiffractorFile, 
						   Position *W1, Position *W2, Position *W3, char **atomLabel, int nAtom)
: Diffractor(DiffractorFile, W1, W2, W3, atomLabel, nAtom)
{
newM(v1, NIntensity, NAtom);
newM(v2, NIntensity, NAtom);
}


Pbar1Diffractor::~Pbar1Diffractor()
{
deleteM(v1);
deleteM(v2);
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////


void Pbar1Diffractor::TransformDerivatives(double *IDerivative, 
										Position *XDerivative)
{
int    n,m;
double x1/*, x2*/, s1/*, s2*/, s3/*, s4*/, s5;//, s6;

for(n=1;n<=NAtom;n++)	
	for(m=1;m<=3;m++)
		FDeriv[n][m] = 0.0;


double *w1;//,*w2;

for(m=1;m<=NIntensity;m++)
	{
	x1 = -IDerivative[m]*ReAmplitude[m];
//	x2 =  IDerivative[m]*ImAmplitude[m];

	s1 = DH[m] * x1;
//	s2 = DH[m] * x2;
	s3 = DK[m] * x1;
//	s4 = DK[m] * x2;
	s5 = DL[m] * x1;
//	s6 = DL[m] * x2;

	w1 = v1[m];
//	w2 = v2[m];

	for(n=1;n<=NAtom;n++)
		{
		FDeriv[n][1] += w1[n]*s1;// + w2[n]*s2; 
		FDeriv[n][2] += w1[n]*s3;// + w2[n]*s4; 
		FDeriv[n][3] += w1[n]*s5;// + w2[n]*s6; 		
		}
	}

MakeXDerivative(XDerivative);	
}


void Pbar1Diffractor::FourierTransform(Position *AtomPosition)
{
int    n,m;
double SinKYplusLZ,CosKYplusLZ;

double s,c;

int h,k,l;

for(n=1;n<=NIntensity;n++)
	{
	ReAmplitude[n] = 0.0;	
	ImAmplitude[n] = 0.0;
	}

for(m=1;m<=NAtom;m++)
	{
	SetUpTrig(AtomPosition[m]);

	for(n=1;n<=NIntensity;n++)
		{
		h = H[n];
		k = K[n];
		l = L[n];

		SinKYplusLZ	= SinY[k]*CosZ[l] + CosY[k]*SinZ[l];
		CosKYplusLZ = CosY[k]*CosZ[l] - SinY[k]*SinZ[l];

		s  = SinX[h]*CosKYplusLZ + CosX[h]*SinKYplusLZ;
		c  = CosX[h]*CosKYplusLZ - SinX[h]*SinKYplusLZ;

		s *= FiniteFactor[m][n];
		c *= FiniteFactor[m][n];

		ReAmplitude[n] += 2.0*c;  
	//	ImAmplitude[n] += s;

		v1[n][m]  = s;
		//v2[n][m]  = c;
		}
	}
}


